Artificial Intelligence untuk Nowcasting

Package

#|output: false
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.4     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readxl)
library(timetk)
library(skimr)
library(torch)

df <- read_excel("Data_PDRB_siap.xlsx")
glimpse(df)

Rows: 108
Columns: 9
$ date              <chr> "2015-01-01", "2015-02-01", "2015-03-01", "2015-04-0…
$ ADHK              <dbl> NA, NA, 20908.25, NA, NA, 21772.33, NA, NA, 24374.33…
$ prod_tbg          <dbl> 55258.07, 60589.93, 65392.61, 65309.37, 60522.22, 67…
$ kredit_konstruksi <dbl> 415006364126, 425296793497, 427231593289, 4268685445…
$ ekspor            <dbl> 63006048, 93951704, 95186751, 2977944, 287826539, 85…
$ pmi_jepang        <dbl> 52.2, 51.6, 50.3, 49.9, 50.9, 50.1, 51.2, 51.7, 51.0…
$ ihk               <dbl> 98.85237, 98.85837, 98.85527, 98.85477, 98.85197, 98…
$ soi               <dbl> -6.8, -7.8, 0.6, -11.3, -3.1, -13.7, -12.0, -14.7, -…
$ ikk               <dbl> 128.80, 123.50, 118.90, 114.30, 118.20, 121.10, 124.…

skim_without_charts(df)

Data summary

Name	df
Number of rows	108
Number of columns	9
_______________________
Column type frequency:
character	1
numeric	8
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
date	0	1	10	10	0	108	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100
ADHK	72	0.33	2.381665e+04	1.471340e+03	2.090825e+04	2.286731e+04	2.353550e+04	2.461094e+04	2.686253e+04
prod_tbg	0	1.00	5.130405e+04	2.225696e+04	1.034700e+04	3.086954e+04	4.929497e+04	6.864869e+04	1.002820e+05
kredit_konstruksi	0	1.00	1.032414e+12	3.340913e+11	4.150064e+11	7.217620e+11	1.146078e+12	1.276832e+12	1.572429e+12
ekspor	0	1.00	1.093068e+08	9.876552e+07	2.098050e+05	3.590831e+07	7.953155e+07	1.575475e+08	4.055064e+08
pmi_jepang	0	1.00	5.069000e+01	2.780000e+00	3.840000e+01	4.907000e+01	5.110000e+01	5.270000e+01	5.540000e+01
ihk	0	1.00	1.035400e+02	5.290000e+00	9.876000e+01	9.883000e+01	1.026100e+02	1.052700e+02	1.159700e+02
soi	2	0.98	1.090000e+00	1.002000e+01	-2.200000e+01	-5.750000e+00	1.450000e+00	8.170000e+00	2.260000e+01
ikk	0	1.00	1.123600e+02	1.353000e+01	7.558000e+01	1.068600e+02	1.141500e+02	1.222100e+02	1.364200e+02

df <- mutate(df,date = ymd(date))
glimpse(df)

Rows: 108
Columns: 9
$ date              <date> 2015-01-01, 2015-02-01, 2015-03-01, 2015-04-01, 201…
$ ADHK              <dbl> NA, NA, 20908.25, NA, NA, 21772.33, NA, NA, 24374.33…
$ prod_tbg          <dbl> 55258.07, 60589.93, 65392.61, 65309.37, 60522.22, 67…
$ kredit_konstruksi <dbl> 415006364126, 425296793497, 427231593289, 4268685445…
$ ekspor            <dbl> 63006048, 93951704, 95186751, 2977944, 287826539, 85…
$ pmi_jepang        <dbl> 52.2, 51.6, 50.3, 49.9, 50.9, 50.1, 51.2, 51.7, 51.0…
$ ihk               <dbl> 98.85237, 98.85837, 98.85527, 98.85477, 98.85197, 98…
$ soi               <dbl> -6.8, -7.8, 0.6, -11.3, -3.1, -13.7, -12.0, -14.7, -…
$ ikk               <dbl> 128.80, 123.50, 118.90, 114.30, 118.20, 121.10, 124.…

Grafik Time Series

plot_time_series(.data=na.omit(df),
                 .date_var = date,
                 .value = ADHK,
                 .interactive = TRUE,
                 .title =  "ADHK",
                 .smooth = FALSE)

fungsi pembantu

multi_plot_time_series <- function(data,date,exclude_var=NULL,.interactive=TRUE,n_col=2,n_row=2,.title="Multiple Time Series"){
data %>% 
  select(-all_of(exclude_var)) %>% 
  select(all_of(date),where(is.numeric)) %>% 
  pivot_longer(cols = -all_of(date),
               names_to = "variable",
               values_to = "value") %>% 
  group_by(variable) %>% 
    plot_time_series(.data=,
                 .date_var = date,
                 .value = value,
                 .interactive = .interactive,
                 .title =  .title,
                 .facet_ncol = n_col,
                 .facet_nrow = n_row,
                 .smooth = FALSE)
}  

multi_plot_time_series(df,
                       date = "date",
                       exclude_var = "ADHK",
                       .interactive=FALSE,
                       n_col = 2)

Grafik Time Series Differencing

tk_augment_differences(.data=na.omit(df),
                        .value = ADHK,
                        .differences = 1) %>% 
  plot_time_series(.date_var = date,
                   .value = ADHK_lag1_diff1,
                   .interactive = TRUE,
                   .title =  "ADHK",
                   .smooth = FALSE)

Multi Input Multi Output (MIMO)

Fungsi pembantu

prepare_mimo_data <- function(data,
                              date_col,
                              input_vars,
                              output_vars,
                              lags       = 1:12,
                              horizon    = 1,
                              remove_na  = TRUE) {
  
  # Tidyeval the date column
  date_col <- rlang::enquo(date_col)
  
  # 1) Order by time index
  df_prep <- data %>%
    dplyr::arrange(!!date_col)
  
  # 2) Generate lagged inputs via timetk
  #    Creates columns like: sales_lag1, sales_lag2, ..., price_lag1, ...
  df_prep <- df_prep %>%
    timetk::tk_augment_lags(
      .value = all_of(input_vars),
      .lags  = lags
    )
  # 3) Generate future targets via dplyr::lead()
  #    Creates columns like: sales_h1, sales_h2, ...
  df_prep <- df_prep %>%
    timetk::tk_augment_leads(
      .value = all_of(output_vars),
      .lags  = -horizon
    )
  
      # Build vector of all generated column names
    lag_cols    <- df_prep %>% select(contains("lag")) %>% names()
    lead_cols   <- df_prep %>% select(contains("lead")) %>% names()
    all_new_cols <- c(sort(lag_cols,decreasing = TRUE), lead_cols)
  # 4) Optionally drop rows with NAs in any of the new columns
  if (remove_na) {
    
    df_prep <- df_prep %>%
      tidyr::drop_na(dplyr::all_of(all_new_cols))
  }
  
  # Return the prepared tibble
  df_prep <- df_prep %>% 
              dplyr::select(!!date_col,
                     dplyr::all_of(all_new_cols)) %>% 
              dplyr::rename("date_lg0"=!!date_col)
  #nm_df_prep <- df_prep %>% select(-!!date_col) %>% names()
  #date_nm <- df_prep %>% select(!!date_col) %>% names()
  #names(df_prep) <- c(date_nm,sort(nm_df_prep,decreasing = FALSE))
  return(df_prep)
}

Struktur Data MIMO

prepare_mimo_data(data = na.omit(df),
                  date_col = date,
                  input_vars = c("ADHK"),
                  output_vars = c("ADHK"),
                  lags = 0:2,
                  remove_na = FALSE,
                  horizon = 1:4)

# A tibble: 35 × 8
   date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
   <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
 1 2015-03-01       NA        NA     20908.     21772.     24374.     22283.
 2 2015-06-01       NA     20908.    21772.     24374.     22283.     22643.
 3 2015-09-01    20908.    21772.    24374.     22283.     22643.     23601.
 4 2015-12-01    21772.    24374.    22283.     22643.     23601.     23122.
 5 2016-03-01    24374.    22283.    22643.     23601.     23122.     21901.
 6 2016-06-01    22283.    22643.    23601.     23122.     21901.     23228.
 7 2016-12-01    22643.    23601.    23122.     21901.     23228.     26217.
 8 2017-03-01    23601.    23122.    21901.     23228.     26217.     23263.
 9 2017-06-01    23122.    21901.    23228.     26217.     23263.     21902.
10 2017-09-01    21901.    23228.    26217.     23263.     21902.     22931.
# ℹ 25 more rows
# ℹ 1 more variable: ADHK_lead4 <dbl>

LSTM

Fungsi pembantu LSTM

train_lstm_mimo <- function(data,
                            date_col,
                            input_cols,
                            output_cols,
                            val_split    = 0.2,
                            epochs       = 50,
                            batch_size   = 32,
                            lr           = 1e-3,
                            optimizer    = c("adam","sgd"),
                            hidden_size  = 50,
                            num_layers   = 1,
                            activation   = c("tanh","relu","linear"),
                            dropout      = 0.0,
                            weight_decay = 0.0) {
  optimizer <- match.arg(optimizer)
  activation <- match.arg(activation)
  date_col   <- rlang::ensym(date_col)

  # 1) Order by time index
  data <- data %>% arrange(!!date_col)

  # 2) Split data
  n     <- nrow(data)
  n_val <- floor(val_split * n)
  train_df <- data[1:(n - n_val), ]
  val_df   <- data[(n - n_val + 1):n, ]

  # 3) Compute robust scaler on train_df
  input_median  <- sapply(input_cols, function(col) median(train_df[[col]], na.rm = TRUE))
  input_iqr     <- sapply(input_cols, function(col) IQR(train_df[[col]], na.rm = TRUE))
  output_median <- sapply(output_cols,function(col) median(train_df[[col]], na.rm = TRUE))
  output_iqr    <- sapply(output_cols,function(col) IQR(train_df[[col]], na.rm = TRUE))
  scaler <- list(
    input_median  = input_median,
    input_iqr     = input_iqr,
    output_median = output_median,
    output_iqr    = output_iqr
  )

  # 4) Apply scaling to train and validation sets
  for (col in input_cols) {
    train_df[[col]] <- (train_df[[col]] - scaler$input_median[col]) / scaler$input_iqr[col]
    val_df[[col]]   <- (val_df[[col]]   - scaler$input_median[col]) / scaler$input_iqr[col]
  }
  for (col in output_cols) {
    train_df[[col]] <- (train_df[[col]] - scaler$output_median[col]) / scaler$output_iqr[col]
    val_df[[col]]   <- (val_df[[col]]   - scaler$output_median[col]) / scaler$output_iqr[col]
  }

  # 5) Define the LSTM module
  LSTMModel <- nn_module(
    "LSTMModel",
    initialize = function(input_size, hidden_size, num_layers, dropout, output_size, activation) {
      self$lstm <- nn_lstm(
        input_size  = input_size,
        hidden_size = hidden_size,
        num_layers  = num_layers,
        batch_first = TRUE,
        dropout     = dropout
      )
      self$fc  <- nn_linear(hidden_size, output_size)
      self$act <- switch(
        activation,
        tanh   = nn_tanh(),
        relu   = nn_relu(),
        linear = nn_identity()
      )
    },
    forward = function(x) {
      out    <- self$lstm(x)
      h_last <- out[[1]][ , dim(out[[1]])[2], ]
      h_act  <- self$act(h_last)
      self$fc(h_act)
    }
  )

  # 6) Prepare torch datasets
  make_ds <- function(df) {
    x_mat <- as.matrix(df[, input_cols])
    y_mat <- as.matrix(df[, output_cols])
    X <- torch_tensor(x_mat, dtype = torch_float())$view(c(nrow(x_mat), -1, length(input_cols)))
    Y <- torch_tensor(y_mat, dtype = torch_float())
    list(x = X, y = Y)
  }
  train_ds <- make_ds(train_df)
  val_ds   <- make_ds(val_df)

  # 7) Instantiate model and optimizer
  model <- LSTMModel(
    input_size  = length(input_cols),
    hidden_size = hidden_size,
    num_layers  = num_layers,
    dropout     = dropout,
    output_size = length(output_cols),
    activation  = activation
  )
  optim <- switch(
    optimizer,
    adam = optim_adam(model$parameters, lr = lr, weight_decay = weight_decay),
    sgd  = optim_sgd(model$parameters, lr = lr, weight_decay = weight_decay)
  )
  criterion <- nn_mse_loss()

  # 8) Training loop
  train_loss <- numeric(epochs)
  val_loss   <- numeric(epochs)
  for (e in seq_len(epochs)) {
    model$train()
    optim$zero_grad()
    preds_train <- model(train_ds$x)
    loss_train  <- criterion(preds_train, train_ds$y)
    loss_train$backward()
    optim$step()
    train_loss[e] <- loss_train$item()

    model$eval()
    with_no_grad({
      preds_val    <- model(val_ds$x)
      val_loss[e]  <- criterion(preds_val, val_ds$y)$item()
    })
  }

  list(
    model       = model,
    train_loss  = train_loss,
    val_loss    = val_loss,
    scaler      = scaler,
    input_cols  = input_cols,
    output_cols = output_cols,
    date_col    = rlang::as_string(date_col)
  )
}

plot_lstm_history <- function(history) {
  df <- tibble(
    epoch = seq_along(history$train_loss),
    training = sqrt(history$train_loss),
    validation = sqrt(history$val_loss)
  ) %>%
    pivot_longer(-epoch, names_to = "data", values_to = "loss")

  ggplot(df, aes(epoch, loss, color = data)) +
    geom_line(size = 1) +
    labs(
      title = "Training vs Validation Loss",
      x     = "Epoch",
      y     = "RMSE"
    ) +
    theme_minimal()
}

predict_lstm <- function(history, new_data) {
  model       <- history$model
  scaler      <- history$scaler
  input_cols  <- history$input_cols
  output_cols <- history$output_cols
  date_col    <- history$date_col

  dates <- new_data[[date_col]]
  x_mat <- as.matrix(new_data[, input_cols])
  X     <- torch_tensor(x_mat, dtype = torch_float())$view(c(nrow(x_mat), -1, ncol(x_mat)))

  model$eval()
  with_no_grad({ pred_scaled <- model(X) })
  pred_scaled_mat <- as.matrix(pred_scaled)

  # Inverse robust scaling
  pred_orig <- sweep(pred_scaled_mat, 2, scaler$output_iqr, `*`)
  pred_orig <- sweep(pred_orig, 2, scaler$output_median, `+`)

  # Build output tibble
  pred_df <- as_tibble(pred_orig)
  names(pred_df) <- output_cols
  tibble(date = dates) %>% bind_cols(pred_df)
}

forecast_lstm <- function(history,output_var,train_data,test_data=NULL){
  forecast <- predict_lstm(history = history,
                           new_data = train_data %>% slice_tail(n = 1)) %>% 
                            pivot_longer(-where(is.Date),
                                         values_to = "forecast"
                                         )

new_dates <-  train_data %>% 
              slice_tail(n = 1) %>% 
              select(where(is.Date)) %>% 
              pull()

for (i in seq(nrow(forecast))) {
  new_dates[i+1] <- ymd(new_dates[i]) + months(3)
}
forecast_df <- tibble(date=new_dates[-1],
                      forecast = forecast %>% pull(forecast)
                      )
result <- train_data %>% 
    select(where(is.Date),all_of(output_var)) %>% 
    pivot_longer(-where(is.Date),names_to = "type") %>% 
    mutate(type="actual") %>% 
    rename_with(.cols = where(is.Date),.fn = function(x) "date" ) %>% 
    bind_rows(
forecast_df %>% 
  pivot_longer(-where(is.Date),names_to = "type")
    ) %>% 
  rename_with(.cols=where(is.Date),.fn =  function(x) "date" )

if(!is.null(test_data)){
test_data <- test_data %>% 
              select(where(is.Date),all_of(output_var)) %>% 
              pivot_longer(-where(is.Date),names_to = "type") %>% 
              mutate(type="actual") %>% 
              rename_with(.cols = where(is.Date),.fn = function(x) "date" )
result <- bind_rows(result,test_data)
}

return(result)
}

Satu Variabel

mimo_df1 <- prepare_mimo_data(data = na.omit(df),
                  date_col = date,
                  input_vars = c("ADHK"),
                  output_vars = c("ADHK"),
                  lags = 0:2,
                  remove_na = TRUE,
                  horizon = 1:4)
mimo_df1

# A tibble: 29 × 8
   date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
   <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
 1 2015-09-01    20908.    21772.    24374.     22283.     22643.     23601.
 2 2015-12-01    21772.    24374.    22283.     22643.     23601.     23122.
 3 2016-03-01    24374.    22283.    22643.     23601.     23122.     21901.
 4 2016-06-01    22283.    22643.    23601.     23122.     21901.     23228.
 5 2016-12-01    22643.    23601.    23122.     21901.     23228.     26217.
 6 2017-03-01    23601.    23122.    21901.     23228.     26217.     23263.
 7 2017-06-01    23122.    21901.    23228.     26217.     23263.     21902.
 8 2017-09-01    21901.    23228.    26217.     23263.     21902.     22931.
 9 2017-12-01    23228.    26217.    23263.     21902.     22931.     22560.
10 2018-03-01    26217.    23263.    21902.     22931.     22560.     22955.
# ℹ 19 more rows
# ℹ 1 more variable: ADHK_lead4 <dbl>

train_df <- mimo_df1 %>% 
            filter(date_lg0<"2022-01-01")
train_df

# A tibble: 25 × 8
   date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
   <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
 1 2015-09-01    20908.    21772.    24374.     22283.     22643.     23601.
 2 2015-12-01    21772.    24374.    22283.     22643.     23601.     23122.
 3 2016-03-01    24374.    22283.    22643.     23601.     23122.     21901.
 4 2016-06-01    22283.    22643.    23601.     23122.     21901.     23228.
 5 2016-12-01    22643.    23601.    23122.     21901.     23228.     26217.
 6 2017-03-01    23601.    23122.    21901.     23228.     26217.     23263.
 7 2017-06-01    23122.    21901.    23228.     26217.     23263.     21902.
 8 2017-09-01    21901.    23228.    26217.     23263.     21902.     22931.
 9 2017-12-01    23228.    26217.    23263.     21902.     22931.     22560.
10 2018-03-01    26217.    23263.    21902.     22931.     22560.     22955.
# ℹ 15 more rows
# ℹ 1 more variable: ADHK_lead4 <dbl>

test_df <- mimo_df1 %>% 
            filter(date_lg0>="2022-01-01")
test_df

# A tibble: 4 × 8
  date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
  <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
1 2022-03-01    24341.    24211.    24428.     25662.     26070.     25914.
2 2022-06-01    24211.    24428.    25662.     26070.     25914.     25296.
3 2022-09-01    24428.    25662.    26070.     25914.     25296.     25267.
4 2022-12-01    25662.    26070.    25914.     25296.     25267.     26480.
# ℹ 1 more variable: ADHK_lead4 <dbl>

input_cols <- names(select(train_df,contains("lag")))
input_cols

[1] "ADHK_lag2" "ADHK_lag1" "ADHK_lag0"

output_cols <- names(select(train_df,contains("lead")))
output_cols

[1] "ADHK_lead1" "ADHK_lead2" "ADHK_lead3" "ADHK_lead4"

mod <- train_lstm_mimo(data = train_df,
                     input_cols = input_cols,
                     output_cols = output_cols,
                     date_col = "date_lg0",
                     val_split = 0.2,
                     epochs = 75,
                     batch_size = 32,
                     optimizer = "adam",
                     hidden_size = 10,
                     num_layers = 1,
                     activation = "linear")
mod$model

An `nn_module` containing 644 parameters.

── Modules ─────────────────────────────────────────────────────────────────────
• lstm: <nn_lstm> #600 parameters
• fc: <nn_linear> #44 parameters
• act: <nn_identity> #0 parameters

plot_lstm_history(mod)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

res1 <- forecast_lstm(history = mod,
                      output_var = "ADHK_lag0",
                      train_data = train_df,
                      test_data = test_df)

Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.

res1 %>% filter(type=="forecast")

# A tibble: 4 × 3
  date       type      value
  <date>     <chr>     <dbl>
1 2022-03-01 forecast 23281.
2 2022-06-01 forecast 22984.
3 2022-09-01 forecast 23260.
4 2022-12-01 forecast 23472.

# RMSE
yardstick::rmse_vec(truth = filter(res1,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res1,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 2363.496

# MAPE
yardstick::mape_vec(truth = filter(res1,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res1,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 8.833074

plot_time_series(.data=res1,
                 .date_var = date,
                 .value = value,
                 .interactive = TRUE,
                 .title =  "Forecast Plot",
                 .color_var = type,
                 .legend_show = FALSE,
                 .smooth = FALSE)

Multi Variabel

#input_vars <- names(select(df,-c(date,ADHK)))
input_vars <- c("ekspor","ihk")
input_vars

[1] "ekspor" "ihk"   

mimo_df2 <- prepare_mimo_data(data = na.omit(df),
                  date_col = date,
                  input_vars = c("ADHK",input_vars),
                  output_vars = c("ADHK"),
                  lags = 0:2,
                  remove_na = TRUE,
                  horizon = 1:4)
mimo_df2

# A tibble: 29 × 14
   date_lg0   ihk_lag2 ihk_lag1 ihk_lag0 ekspor_lag2 ekspor_lag1 ekspor_lag0
   <date>        <dbl>    <dbl>    <dbl>       <dbl>       <dbl>       <dbl>
 1 2015-09-01     98.9     98.9     98.8    95186751    85868149   146184552
 2 2015-12-01     98.9     98.8     98.8    85868149   146184552   277363045
 3 2016-03-01     98.8     98.8     98.8   146184552   277363045    80837466
 4 2016-06-01     98.8     98.8     98.8   277363045    80837466   256450298
 5 2016-12-01     98.8     98.8     98.8    80837466   256450298   195677273
 6 2017-03-01     98.8     98.8     98.8   256450298   195677273   135640406
 7 2017-06-01     98.8     98.8     98.8   195677273   135640406    59935547
 8 2017-09-01     98.8     98.8     98.8   135640406    59935547    75541458
 9 2017-12-01     98.8     98.8     98.8    59935547    75541458   158378320
10 2018-03-01     98.8     98.8     98.9    75541458   158378320    43910678
# ℹ 19 more rows
# ℹ 7 more variables: ADHK_lag2 <dbl>, ADHK_lag1 <dbl>, ADHK_lag0 <dbl>,
#   ADHK_lead1 <dbl>, ADHK_lead2 <dbl>, ADHK_lead3 <dbl>, ADHK_lead4 <dbl>

train_df <- mimo_df2 %>% 
            filter(date_lg0<"2022-01-01")
train_df

# A tibble: 25 × 14
   date_lg0   ihk_lag2 ihk_lag1 ihk_lag0 ekspor_lag2 ekspor_lag1 ekspor_lag0
   <date>        <dbl>    <dbl>    <dbl>       <dbl>       <dbl>       <dbl>
 1 2015-09-01     98.9     98.9     98.8    95186751    85868149   146184552
 2 2015-12-01     98.9     98.8     98.8    85868149   146184552   277363045
 3 2016-03-01     98.8     98.8     98.8   146184552   277363045    80837466
 4 2016-06-01     98.8     98.8     98.8   277363045    80837466   256450298
 5 2016-12-01     98.8     98.8     98.8    80837466   256450298   195677273
 6 2017-03-01     98.8     98.8     98.8   256450298   195677273   135640406
 7 2017-06-01     98.8     98.8     98.8   195677273   135640406    59935547
 8 2017-09-01     98.8     98.8     98.8   135640406    59935547    75541458
 9 2017-12-01     98.8     98.8     98.8    59935547    75541458   158378320
10 2018-03-01     98.8     98.8     98.9    75541458   158378320    43910678
# ℹ 15 more rows
# ℹ 7 more variables: ADHK_lag2 <dbl>, ADHK_lag1 <dbl>, ADHK_lag0 <dbl>,
#   ADHK_lead1 <dbl>, ADHK_lead2 <dbl>, ADHK_lead3 <dbl>, ADHK_lead4 <dbl>

test_df <- mimo_df2 %>% 
            filter(date_lg0>="2022-01-01")
test_df

# A tibble: 4 × 14
  date_lg0   ihk_lag2 ihk_lag1 ihk_lag0 ekspor_lag2 ekspor_lag1 ekspor_lag0
  <date>        <dbl>    <dbl>    <dbl>       <dbl>       <dbl>       <dbl>
1 2022-03-01     105.     106.     108.   101489132    58472552   212100511
2 2022-06-01     106.     108.     111.    58472552   212100511   323431909
3 2022-09-01     108.     111.     112.   212100511   323431909   212074483
4 2022-12-01     111.     112.     113.   323431909   212074483   232322051
# ℹ 7 more variables: ADHK_lag2 <dbl>, ADHK_lag1 <dbl>, ADHK_lag0 <dbl>,
#   ADHK_lead1 <dbl>, ADHK_lead2 <dbl>, ADHK_lead3 <dbl>, ADHK_lead4 <dbl>

input_cols <- names(select(train_df,contains("lag")))
input_cols

[1] "ihk_lag2"    "ihk_lag1"    "ihk_lag0"    "ekspor_lag2" "ekspor_lag1"
[6] "ekspor_lag0" "ADHK_lag2"   "ADHK_lag1"   "ADHK_lag0"  

output_cols <- names(select(train_df,contains("lead")))
output_cols

[1] "ADHK_lead1" "ADHK_lead2" "ADHK_lead3" "ADHK_lead4"

mod2 <- train_lstm_mimo(data = train_df,
                     input_cols = input_cols,
                     output_cols = output_cols,
                     date_col = "date_lg0",
                     val_split = 0.2,
                     epochs = 100,
                     batch_size = 32,
                     optimizer = "adam",
                     hidden_size = 100,
                     num_layers = 1,
                     activation = "linear")
mod2$model

An `nn_module` containing 44,804 parameters.

── Modules ─────────────────────────────────────────────────────────────────────
• lstm: <nn_lstm> #44,400 parameters
• fc: <nn_linear> #404 parameters
• act: <nn_identity> #0 parameters

plot_lstm_history(mod2)

res2 <- forecast_lstm(history = mod2,
                      output_var = "ADHK_lag0",
                      train_data = train_df,
                      test_data = test_df)
res2 %>% filter(type=="forecast")

# A tibble: 4 × 3
  date       type      value
  <date>     <chr>     <dbl>
1 2022-03-01 forecast 24462.
2 2022-06-01 forecast 26022.
3 2022-09-01 forecast 22610.
4 2022-12-01 forecast 24235.

# RMSE
yardstick::rmse_vec(truth = filter(res2,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res2,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 1931.565

# MAPE
yardstick::mape_vec(truth = filter(res2,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res2,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 5.324029

plot_time_series(.data=res2,
                 .date_var = date,
                 .value = value,
                 .interactive = TRUE,
                 .title =  "Forecast Plot",
                 .color_var = type,
                 .legend_show = FALSE,
                 .smooth = FALSE)

RNN

Fungsi Pembantu

train_rnn_mimo <- function(data,
                            date_col,
                            input_cols,
                            output_cols,
                            val_split    = 0.2,
                            epochs       = 50,
                            batch_size   = 32,
                            lr           = 1e-3,
                            optimizer    = c("adam","sgd"),
                            hidden_size  = 50,
                            num_layers   = 1,
                            activation   = c("tanh","relu","linear"),
                            dropout      = 0.0,
                            weight_decay = 0.0) {
  optimizer <- match.arg(optimizer)
  activation <- match.arg(activation)
  date_col   <- rlang::ensym(date_col)

  # 1) Order by time index
  data <- data %>% arrange(!!date_col)

  # 2) Split data
  n     <- nrow(data)
  n_val <- floor(val_split * n)
  train_df <- data[1:(n - n_val), ]
  val_df   <- data[(n - n_val + 1):n, ]

  # 3) Compute robust scaler on train_df
  input_median  <- sapply(input_cols, function(col) median(train_df[[col]], na.rm = TRUE))
  input_iqr     <- sapply(input_cols, function(col) IQR(train_df[[col]], na.rm = TRUE))
  output_median <- sapply(output_cols,function(col) median(train_df[[col]], na.rm = TRUE))
  output_iqr    <- sapply(output_cols,function(col) IQR(train_df[[col]], na.rm = TRUE))
  scaler <- list(
    input_median  = input_median,
    input_iqr     = input_iqr,
    output_median = output_median,
    output_iqr    = output_iqr
  )

  # 4) Apply scaling to train and validation sets
  for (col in input_cols) {
    train_df[[col]] <- (train_df[[col]] - scaler$input_median[col]) / scaler$input_iqr[col]
    val_df[[col]]   <- (val_df[[col]]   - scaler$input_median[col]) / scaler$input_iqr[col]
  }
  for (col in output_cols) {
    train_df[[col]] <- (train_df[[col]] - scaler$output_median[col]) / scaler$output_iqr[col]
    val_df[[col]]   <- (val_df[[col]]   - scaler$output_median[col]) / scaler$output_iqr[col]
  }

  # 5) Define the RNN module
  RNNModel <- nn_module(
  "RNNModel",
  initialize = function(input_size,
                        hidden_size,
                        num_layers,
                        dropout,
                        output_size,
                        activation,
                        nonlinearity = c("tanh","relu")) {
    nonlinearity <- match.arg(nonlinearity)
    self$rnn <- nn_rnn(
      input_size   = input_size,
      hidden_size  = hidden_size,
      num_layers   = num_layers,
      nonlinearity = nonlinearity,
      batch_first  = TRUE,
      dropout      = dropout
    )
    self$fc  <- nn_linear(hidden_size, output_size)
    self$act <- switch(
      activation,
      tanh   = nn_tanh(),
      relu   = nn_relu(),
      linear = nn_identity()
    )
  },
  forward = function(x) {
    # x: [batch, seq_len, input_size]
    out   <- self$rnn(x)
    # out[[1]] is the output at every time step: shape [batch, seq_len, hidden]
    last  <- out[[1]][ , dim(out[[1]])[2], ]
    h_act <- self$act(last)
    self$fc(h_act)
  }
)


  # 6) Prepare torch datasets
  make_ds <- function(df) {
    x_mat <- as.matrix(df[, input_cols])
    y_mat <- as.matrix(df[, output_cols])
    X <- torch_tensor(x_mat, dtype = torch_float())$view(c(nrow(x_mat), -1, length(input_cols)))
    Y <- torch_tensor(y_mat, dtype = torch_float())
    list(x = X, y = Y)
  }
  train_ds <- make_ds(train_df)
  val_ds   <- make_ds(val_df)

  # 7) Instantiate model and optimizer
  model <- RNNModel(
  input_size   = length(input_cols),
  hidden_size  = hidden_size,
  num_layers   = num_layers,
  dropout      = dropout,
  output_size  = length(output_cols),
  activation   = activation,
  nonlinearity = "relu")

  optim <- switch(
    optimizer,
    adam = optim_adam(model$parameters, lr = lr, weight_decay = weight_decay),
    sgd  = optim_sgd(model$parameters, lr = lr, weight_decay = weight_decay)
  )
  criterion <- nn_mse_loss()

  # 8) Training loop
  train_loss <- numeric(epochs)
  val_loss   <- numeric(epochs)
  for (e in seq_len(epochs)) {
    model$train()
    optim$zero_grad()
    preds_train <- model(train_ds$x)
    loss_train  <- criterion(preds_train, train_ds$y)
    loss_train$backward()
    optim$step()
    train_loss[e] <- loss_train$item()

    model$eval()
    with_no_grad({
      preds_val    <- model(val_ds$x)
      val_loss[e]  <- criterion(preds_val, val_ds$y)$item()
    })
  }

  list(
    model       = model,
    train_loss  = train_loss,
    val_loss    = val_loss,
    scaler      = scaler,
    input_cols  = input_cols,
    output_cols = output_cols,
    date_col    = rlang::as_string(date_col)
  )
}

plot_rnn_history <- function(history) {
  df <- tibble(
    epoch = seq_along(history$train_loss),
    training = sqrt(history$train_loss),
    validation = sqrt(history$val_loss)
  ) %>%
    pivot_longer(-epoch, names_to = "data", values_to = "loss")

  ggplot(df, aes(epoch, loss, color = data)) +
    geom_line(size = 1) +
    labs(
      title = "Training vs Validation Loss",
      x     = "Epoch",
      y     = "RMSE"
    ) +
    theme_minimal()
}

predict_rnn <- function(history, new_data) {
  model       <- history$model
  scaler      <- history$scaler
  input_cols  <- history$input_cols
  output_cols <- history$output_cols
  date_col    <- history$date_col

  dates <- new_data[[date_col]]
  x_mat <- as.matrix(new_data[, input_cols])
  X     <- torch_tensor(x_mat, dtype = torch_float())$view(c(nrow(x_mat), -1, ncol(x_mat)))

  model$eval()
  with_no_grad({ pred_scaled <- model(X) })
  pred_scaled_mat <- as.matrix(pred_scaled)

  # Inverse robust scaling
  pred_orig <- sweep(pred_scaled_mat, 2, scaler$output_iqr, `*`)
  pred_orig <- sweep(pred_orig, 2, scaler$output_median, `+`)

  # Build output tibble
  pred_df <- as_tibble(pred_orig)
  names(pred_df) <- output_cols
  tibble(date = dates) %>% bind_cols(pred_df)
}

forecast_rnn <- function(history,output_var,train_data,test_data=NULL){
  forecast <- predict_rnn(history = history,
                           new_data = train_data %>% slice_tail(n = 1)) %>% 
                            pivot_longer(-where(is.Date),
                                         values_to = "forecast"
                                         )

new_dates <-  train_data %>% 
              slice_tail(n = 1) %>% 
              select(where(is.Date)) %>% 
              pull()

for (i in seq(nrow(forecast))) {
  new_dates[i+1] <- ymd(new_dates[i]) + months(3)
}
forecast_df <- tibble(date=new_dates[-1],
                      forecast = forecast %>% pull(forecast)
                      )
result <- train_data %>% 
    select(where(is.Date),all_of(output_var)) %>% 
    pivot_longer(-where(is.Date),names_to = "type") %>% 
    mutate(type="actual") %>% 
    rename_with(.cols = where(is.Date),.fn = function(x) "date" ) %>% 
    bind_rows(
forecast_df %>% 
  pivot_longer(-where(is.Date),names_to = "type")
    ) %>% 
  rename_with(.cols=where(is.Date),.fn =  function(x) "date" )

if(!is.null(test_data)){
test_data <- test_data %>% 
              select(where(is.Date),all_of(output_var)) %>% 
              pivot_longer(-where(is.Date),names_to = "type") %>% 
              mutate(type="actual") %>% 
              rename_with(.cols = where(is.Date),.fn = function(x) "date" )
result <- bind_rows(result,test_data)
}

return(result)
}

Satu Variabel

mimo_df1 <- prepare_mimo_data(data = na.omit(df),
                  date_col = date,
                  input_vars = c("ADHK"),
                  output_vars = c("ADHK"),
                  lags = 0:2,
                  remove_na = TRUE,
                  horizon = 1:4)
mimo_df1

# A tibble: 29 × 8
   date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
   <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
 1 2015-09-01    20908.    21772.    24374.     22283.     22643.     23601.
 2 2015-12-01    21772.    24374.    22283.     22643.     23601.     23122.
 3 2016-03-01    24374.    22283.    22643.     23601.     23122.     21901.
 4 2016-06-01    22283.    22643.    23601.     23122.     21901.     23228.
 5 2016-12-01    22643.    23601.    23122.     21901.     23228.     26217.
 6 2017-03-01    23601.    23122.    21901.     23228.     26217.     23263.
 7 2017-06-01    23122.    21901.    23228.     26217.     23263.     21902.
 8 2017-09-01    21901.    23228.    26217.     23263.     21902.     22931.
 9 2017-12-01    23228.    26217.    23263.     21902.     22931.     22560.
10 2018-03-01    26217.    23263.    21902.     22931.     22560.     22955.
# ℹ 19 more rows
# ℹ 1 more variable: ADHK_lead4 <dbl>

train_df <- mimo_df1 %>% 
            filter(date_lg0<"2022-01-01")
train_df

# A tibble: 25 × 8
   date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
   <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
 1 2015-09-01    20908.    21772.    24374.     22283.     22643.     23601.
 2 2015-12-01    21772.    24374.    22283.     22643.     23601.     23122.
 3 2016-03-01    24374.    22283.    22643.     23601.     23122.     21901.
 4 2016-06-01    22283.    22643.    23601.     23122.     21901.     23228.
 5 2016-12-01    22643.    23601.    23122.     21901.     23228.     26217.
 6 2017-03-01    23601.    23122.    21901.     23228.     26217.     23263.
 7 2017-06-01    23122.    21901.    23228.     26217.     23263.     21902.
 8 2017-09-01    21901.    23228.    26217.     23263.     21902.     22931.
 9 2017-12-01    23228.    26217.    23263.     21902.     22931.     22560.
10 2018-03-01    26217.    23263.    21902.     22931.     22560.     22955.
# ℹ 15 more rows
# ℹ 1 more variable: ADHK_lead4 <dbl>

test_df <- mimo_df1 %>% 
            filter(date_lg0>="2022-01-01")
test_df

# A tibble: 4 × 8
  date_lg0   ADHK_lag2 ADHK_lag1 ADHK_lag0 ADHK_lead1 ADHK_lead2 ADHK_lead3
  <date>         <dbl>     <dbl>     <dbl>      <dbl>      <dbl>      <dbl>
1 2022-03-01    24341.    24211.    24428.     25662.     26070.     25914.
2 2022-06-01    24211.    24428.    25662.     26070.     25914.     25296.
3 2022-09-01    24428.    25662.    26070.     25914.     25296.     25267.
4 2022-12-01    25662.    26070.    25914.     25296.     25267.     26480.
# ℹ 1 more variable: ADHK_lead4 <dbl>

input_cols <- names(select(train_df,contains("lag")))
input_cols

[1] "ADHK_lag2" "ADHK_lag1" "ADHK_lag0"

output_cols <- names(select(train_df,contains("lead")))
output_cols

[1] "ADHK_lead1" "ADHK_lead2" "ADHK_lead3" "ADHK_lead4"

mod3 <- train_rnn_mimo(data = train_df,
                     input_cols = input_cols,
                     output_cols = output_cols,
                     date_col = "date_lg0",
                     val_split = 0.2,
                     epochs = 75,
                     batch_size = 32,
                     optimizer = "adam",
                     hidden_size = 10,
                     num_layers = 1,
                     activation = "linear")
mod3$model

An `nn_module` containing 194 parameters.

── Modules ─────────────────────────────────────────────────────────────────────
• rnn: <nn_rnn> #150 parameters
• fc: <nn_linear> #44 parameters
• act: <nn_identity> #0 parameters

plot_rnn_history(mod3)

res3 <- forecast_rnn(history = mod3,
                      output_var = "ADHK_lag0",
                      train_data = train_df,
                      test_data = test_df)
res3 %>% filter(type=="forecast")

# A tibble: 4 × 3
  date       type         value
  <date>     <chr>        <dbl>
1 2022-03-01 forecast   780486.
2 2022-06-01 forecast -1305499.
3 2022-09-01 forecast   986766.
4 2022-12-01 forecast -3011343.

# RMSE
yardstick::rmse_vec(truth = filter(res3,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res3,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 1767165

# MAPE
yardstick::mape_vec(truth = filter(res3,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res3,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 5921.959

plot_time_series(.data=res3,
                 .date_var = date,
                 .value = value,
                 .interactive = TRUE,
                 .title =  "Forecast Plot",
                 .color_var = type,
                 .legend_show = FALSE,
                 .smooth = FALSE)

Multi Variabel

#input_vars <- names(select(df,-c(date,ADHK)))
input_vars <- c("ekspor","ihk")
input_vars

[1] "ekspor" "ihk"   

mimo_df2 <- prepare_mimo_data(data = na.omit(df),
                  date_col = date,
                  input_vars = c("ADHK",input_vars),
                  output_vars = c("ADHK"),
                  lags = 0:2,
                  remove_na = TRUE,
                  horizon = 1:4)
mimo_df2

# A tibble: 29 × 14
   date_lg0   ihk_lag2 ihk_lag1 ihk_lag0 ekspor_lag2 ekspor_lag1 ekspor_lag0
   <date>        <dbl>    <dbl>    <dbl>       <dbl>       <dbl>       <dbl>
 1 2015-09-01     98.9     98.9     98.8    95186751    85868149   146184552
 2 2015-12-01     98.9     98.8     98.8    85868149   146184552   277363045
 3 2016-03-01     98.8     98.8     98.8   146184552   277363045    80837466
 4 2016-06-01     98.8     98.8     98.8   277363045    80837466   256450298
 5 2016-12-01     98.8     98.8     98.8    80837466   256450298   195677273
 6 2017-03-01     98.8     98.8     98.8   256450298   195677273   135640406
 7 2017-06-01     98.8     98.8     98.8   195677273   135640406    59935547
 8 2017-09-01     98.8     98.8     98.8   135640406    59935547    75541458
 9 2017-12-01     98.8     98.8     98.8    59935547    75541458   158378320
10 2018-03-01     98.8     98.8     98.9    75541458   158378320    43910678
# ℹ 19 more rows
# ℹ 7 more variables: ADHK_lag2 <dbl>, ADHK_lag1 <dbl>, ADHK_lag0 <dbl>,
#   ADHK_lead1 <dbl>, ADHK_lead2 <dbl>, ADHK_lead3 <dbl>, ADHK_lead4 <dbl>

train_df <- mimo_df2 %>% 
            filter(date_lg0<"2022-01-01")
train_df

# A tibble: 25 × 14
   date_lg0   ihk_lag2 ihk_lag1 ihk_lag0 ekspor_lag2 ekspor_lag1 ekspor_lag0
   <date>        <dbl>    <dbl>    <dbl>       <dbl>       <dbl>       <dbl>
 1 2015-09-01     98.9     98.9     98.8    95186751    85868149   146184552
 2 2015-12-01     98.9     98.8     98.8    85868149   146184552   277363045
 3 2016-03-01     98.8     98.8     98.8   146184552   277363045    80837466
 4 2016-06-01     98.8     98.8     98.8   277363045    80837466   256450298
 5 2016-12-01     98.8     98.8     98.8    80837466   256450298   195677273
 6 2017-03-01     98.8     98.8     98.8   256450298   195677273   135640406
 7 2017-06-01     98.8     98.8     98.8   195677273   135640406    59935547
 8 2017-09-01     98.8     98.8     98.8   135640406    59935547    75541458
 9 2017-12-01     98.8     98.8     98.8    59935547    75541458   158378320
10 2018-03-01     98.8     98.8     98.9    75541458   158378320    43910678
# ℹ 15 more rows
# ℹ 7 more variables: ADHK_lag2 <dbl>, ADHK_lag1 <dbl>, ADHK_lag0 <dbl>,
#   ADHK_lead1 <dbl>, ADHK_lead2 <dbl>, ADHK_lead3 <dbl>, ADHK_lead4 <dbl>

test_df <- mimo_df2 %>% 
            filter(date_lg0>="2022-01-01")
test_df

# A tibble: 4 × 14
  date_lg0   ihk_lag2 ihk_lag1 ihk_lag0 ekspor_lag2 ekspor_lag1 ekspor_lag0
  <date>        <dbl>    <dbl>    <dbl>       <dbl>       <dbl>       <dbl>
1 2022-03-01     105.     106.     108.   101489132    58472552   212100511
2 2022-06-01     106.     108.     111.    58472552   212100511   323431909
3 2022-09-01     108.     111.     112.   212100511   323431909   212074483
4 2022-12-01     111.     112.     113.   323431909   212074483   232322051
# ℹ 7 more variables: ADHK_lag2 <dbl>, ADHK_lag1 <dbl>, ADHK_lag0 <dbl>,
#   ADHK_lead1 <dbl>, ADHK_lead2 <dbl>, ADHK_lead3 <dbl>, ADHK_lead4 <dbl>

input_cols <- names(select(train_df,contains("lag")))
input_cols

[1] "ihk_lag2"    "ihk_lag1"    "ihk_lag0"    "ekspor_lag2" "ekspor_lag1"
[6] "ekspor_lag0" "ADHK_lag2"   "ADHK_lag1"   "ADHK_lag0"  

output_cols <- names(select(train_df,contains("lead")))
output_cols

[1] "ADHK_lead1" "ADHK_lead2" "ADHK_lead3" "ADHK_lead4"

mod4 <- train_rnn_mimo(data = train_df,
                     input_cols = input_cols,
                     output_cols = output_cols,
                     date_col = "date_lg0",
                     val_split = 0.2,
                     epochs = 100,
                     batch_size = 32,
                     optimizer = "adam",
                     hidden_size = 100,
                     num_layers = 1,
                     activation = "linear")
mod4$model

An `nn_module` containing 11,504 parameters.

── Modules ─────────────────────────────────────────────────────────────────────
• rnn: <nn_rnn> #11,100 parameters
• fc: <nn_linear> #404 parameters
• act: <nn_identity> #0 parameters

plot_rnn_history(mod4)

res4 <- forecast_lstm(history = mod4,
                      output_var = "ADHK_lag0",
                      train_data = train_df,
                      test_data = test_df)
res4 %>% filter(type=="forecast")

# A tibble: 4 × 3
  date       type            value
  <date>     <chr>           <dbl>
1 2022-03-01 forecast 23144651920.
2 2022-06-01 forecast 41473219858.
3 2022-09-01 forecast 11820379034.
4 2022-12-01 forecast 11120872198.

# RMSE
yardstick::rmse_vec(truth = filter(res4,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res4,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 25095286334

# MAPE
yardstick::mape_vec(truth = filter(res4,
                                   type=="actual",
                                   date>="2022-01-01") %>% pull(value),
                    estimate = filter(res4,
                                   type=="forecast") %>% pull(value)
                                   )

[1] 86153954

plot_time_series(.data=res4,
                 .date_var = date,
                 .value = value,
                 .interactive = TRUE,
                 .title =  "Forecast Plot",
                 .color_var = type,
                 .legend_show = FALSE,
                 .smooth = FALSE)